#### make formula for felm ####
f_felm <- function(yvar, xvars = 0, fevars = 0, fit = 0, ivars = 0, clustervars = 0) {
  formula <- as.formula(
    paste(
      yvar, "~", 
      paste(xvars, collapse = "+"), "|", 
      paste(fevars, collapse = "+"), "|",
      ifelse(fit == 0, 0, paste("(", fit, "~", paste(ivars, collapse = "+"), ")")), "|",
      paste(clustervars, collapse = "+")  
    )
  )
  return(formula)
}

#### Conley standard errors ####
conley_concern <- function(x, data = eb_climate, conley_threads = 4) {
  conley_se <- ConleySEs(
    reg = x,
    unit = "unit", time = "t",
    lat = "lat", lon = "lon",
    kernel = "bartlett", #"uniform"
    dist_fn = "Haversine", #great circle distance
    dist_cutoff = env_spc, 
    lag_cutoff = env_tc,
    data = data, 
    balanced_pnl = FALSE,
    cores = conley_threads
  )
  x$vcv <- x$robustvcv <- x$clustervcv <- conley_se$Spatial_HAC
  ct <- lmtest::coeftest(x, vcov. = conley_se$Spatial_HAC)
  x$cse <- x$se <- ct[,2] # standard errors
  x$ctval <- x$tval <- ct[,3] # T statistic
  x$cpval <- x$pval <- ct[,4] # p value
  return(x)
}

conley_vote <- function(x, data = vote_climate, conley_threads = 4) {
  conley_se <- ConleySEs(
    reg = x,
    unit = "unit", time = "year",
    lat = "lat", lon = "lon",
    kernel = "bartlett", 
    dist_fn = "Haversine", #great circle distance
    dist_cutoff = gv_spc, 
    lag_cutoff = gv_tc, 
    data = data, 
    balanced_pnl = FALSE,
    cores = conley_threads
  )
  x$vcv <- x$robustvcv <- x$clustervcv <- conley_se$Spatial_HAC
  ct <- lmtest::coeftest(x, vcov. = conley_se$Spatial_HAC)
  x$cse <- x$se <- ct[,2] # standard errors
  x$ctval <- x$tval <- ct[,3] # T statistic
  x$cpval <- x$pval <- ct[,4] # p value
  return(x)
}

#### standardise felm ####
standardise <- 
  function(
    x, data, 
    res_y = TRUE, res_x = TRUE, # residualise y and/or x?
    stan_y = TRUE, stan_x = TRUE, none = FALSE, # standardise?
    stage1 = FALSE # standardise first stage of two stage design?
  ){
    stopifnot( # input checks
      (class(x) == "felm" | stage1 == TRUE), 
      is.logical(res_y), is.logical(res_x), is.logical(none)
    )
    if (stage1 == TRUE) {
      x <- x$stage1
    }
    if (none == TRUE) {return(x)} # dont standardise
    if (res_x == TRUE) { # use data after fixed effects transformation
      sd_x <- apply(x$cX, 2, function(x) {sd(x, na.rm = TRUE)})
    } else { # or before FE
      sd_x <- apply(x$X, 2, function(x) {sd(x, na.rm = TRUE)})
    }
    if (res_y == TRUE){
      sd_y <- sd(x$cY, na.rm = TRUE) 
    } else {
      sd_y <- sd(x$response, na.rm = TRUE)
    }
    if (stan_y == FALSE) {
      sd_y <- 1
    }
    x$coefficients <- x$coefficients * (sd_x / sd_y)
    x$beta <- x$beta * (sd_x / sd_y)
    x$cse <- x$cse * (sd_x / sd_y)
    return(x)
  }

#### export html table ####
quietly <- function(x) { 
  sink(tempfile()) 
  on.exit(sink()) 
  invisible(force(x)) 
}
export_html <- function(
  reg, out,
  cov_labels = names(xv),
  dep_label_include = TRUE,
  dep_labels = c("Environmental concern", "Green vote share"),
  col_labels = NULL,
  col_sep = rep(length(xv), 2),
  omit_stats = c("all"),
  add_lines = list(
    c("Unit fixed effects", rep("×", length(xv)*2)),
    c("Period fixed effects", rep("×", length(xv)*2)),
    c("Season fixed effects", rep("×", length(xv))),
    c("Spatial cutoff (km)", rep(env_spc, length(xv)), rep(gv_spc, length(xv))),
    c("Temporal cutoff (years)", rep(env_tc/12, length(xv)), rep(gv_tc, length(xv))),
    c("Observations", sapply(reg, get_n)),
    c("Overall R²", sapply(reg, get_r2_overall)),
    c("Within R²", sapply(reg, get_r2_within))
  )
){
  quietly(
    stargazer(
      reg,
      out = out,
      type = "html",
      omit.stat = omit_stats,
      dep.var.labels = dep_labels,
      dep.var.labels.include = dep_label_include,
      column.labels = col_labels,
      column.separate = col_sep,
      covariate.labels = cov_labels,
      add.lines = add_lines
    )
  )
}

#### get values from felm objects ####
get_n <- function(x){
  tmp <- x
  return(tmp$N)
}
get_r2_overall <- function(x){
  tmp <- summary(x)
  return(round(tmp$r2, 3))
}
get_r2_within <- function(x, stage1 = FALSE){
  if (stage1 == TRUE) {
    tmp <- summary(x$stage1)
  } else {
    tmp <- summary(x)
  }
  return(round(tmp$P.r.squared, 3))
}

get_f <- function(x){
  tmp <- summary(x)
  return(round(tmp$P.fstat["F"], 3))
}
get_pf <- function(x){
  tmp <- summary(x)
  return(round(tmp$P.fstat["p.F"], 3))
  }
get_f_iv <- function(x){
  return(round(x$stage1$iv1fstat$env_concern_12m["F"], 1))
}
get_pf_iv <- function(x){
  return(round(x$stage1$iv1fstat$env_concern_12m["p.F"], 1))
}

